Temporally disordered granular flow: A model of landslides 
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We propose and study numerically a stochastic cellular automaton model for the dynamics of gran- 
ular materials with temporal disorder representing random variation of the diffusion probability 
1 — around threshold value 1 — no during the course of an avalanche. Combined with the slope 
threshold dynamics, the temporal disorder yields a series of secondary instabilities, resembling those 
in realistic granular slides. When the parameter fio is lower than the critical value ~ 0.4, the 
dynamics is dominated by occasional huge sandslides. For the range of values < /^o < 1 the 
critical steady states occur, which are characterized by multifractal scaling properties of the slide 
distributions and continuously varying critical exponents Tx{fJ.o)- The mass distribution exponent 
for no ~ 0.45 is in agreement with the reported value that characterizes Himalayan sandslides. At 
Mo = A'o the exponents governing distributions of large relaxation events reach numerical values 
which are close to those of parity-conserving universality class, whereas for small avalanches they 
are close to the mean- field exponents. 
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I. INTRODUCTION 

Understanding flow in realistic granular materials ap- 
pears to be an important problem from both practical 
and theoretical point of view [Q, Renewed theoret- 
ical interest in this field has concentrated on the origin 
of scaling that characterizes phenomena in slowly driven 
granular materials: Distributions of avalanches in real- 
istic granular piles ^-|^, stratification ||^, compactifica- 
tion ||l^, etc. The central question is: Do granular piles 
self-organize into critical steady states Q and if so, un- 
der what conditions? Another interesting phenomenon 
related to dynamics of granular materials in nature is 
the landscape evolution due to overland and channel 
flow, which results in fractal topography. The underly- 
ing mechanisms of erosion with spatially and temporally 
varying erosion rates are the subject of intensive discus- 
sion in the literature |]lT| . 

It has been understood that realistic flow in slowly 
driven granular piles depends on many parameters, such 
as shapes and sizes (and masses) of individual beans, 
roughness of contact surfaces, their wetting properties, 
etc. Random (or controlled) variations in some of 
these parameters lead to fluctuations of contact angles 
and force distribution fl^ , nonlinear friction, stochas- 
tic character of diffusion, velocity and convection direc- 
tions, and fluctuations in angle of repose. Unidirectional 
flow — reflecting dependence on gravity — is common in all 
granular materials, as well as occurrence of secondary 
avalanches following the initial instability. Molecular 
dynamic (MD) simulations [ p^ and various cellular au- 
tomata models with stochastic relaxation rules 14 -T^] 
have been useful in describing certain aspects of realis- 
tic granular flow. However, comparison with measured 
avalanche properties has been only qualitative. 

In experiments the most often measured quantity is 



the outflow current J, which is defined as the number of 
grains that leave the system when an avalanche hits its 
lower boundary. The probability distribution of outflow 
current P{J) in the steady state obeys the scaling form 
P( J, L) = L-f^giJL-") with l3 = 2iy when the linear size 
L of the pile support is varied, as found in Ref. Q for 
sandpiles of relatively small sizes. Using silicon dioxide 
sand Rosendahl et al. concluded that small and large 
avalanches behave differently and the distribution P(J) 
shows no simple finite-size scaling. Moreover, avalanche 
statistics was found to vary with the size of grains used. 
Measuring the internal avalanches Bretz et al. have 
also observed that two types od statistics are governing 
small and large avalanches. The measured distribution of 
avalanche size exhibits a power-law behavior D{s) ~ 5""^° 
with 1^ « 2.14, which probably applies for avalanches 
of small sizes. The two time (and size) scales were more 
clearly demonstrated recently by MD simulations [p^ , 
leading to two exponents Tg = 2 for short, and Tg = 1.5 
for long time scale. A sophisticated measurement of the 
internal avalanches was done with a one-dimensional ri- 
cepile [Q, in which elongated rice grains were used to 
suppress inertial effects. Scaling properties of the distri- 
bution of dissipated energy were determined, indicating 
that details of the dissipation are responsible for the oc- 
currence of the critical state. In another experiment the 
transport of individual grains was monitored, and the dis- 
tribution of transit time was also found to exhibit robust 
scaling behavior ^ . 

The collected data for the landslides in nature, trig- 
gered by various mechanisms, also exhibit a power-law 
behavior The exponents for the area of slides have 
been estimated in the range Tg = 1.16 — 2.25 ||l^, depend- 
ing on the dominating triggering mechanism and region 
where the data were collected. The distribution of the 
mass collected from Himalayan sandslides is character- 
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ized by the exponent = 0.19 — 0.23 | |T^ . 

In the present work we introduce a new stochastic 
model of directional flow on the two-dimensional square 
lattice in which numerous after-avalanches are generated 
within a certain correlation time due to temporal disor- 
der in the diffusion term. The dynamic rules are a combi- 
nation of stochastic diffusion and deterministic branch- 
ing processes. The diffusion probabilities change ran- 
domly in time, but are space-independent. Fluctuations 
in diffusion probability 1 — ^{t) around threshold value 

I — /zo, which depends on external conditions and thus 
appears as a control parameter, is motivated by fluctua- 
tions in wetting and drying conditions after an avalanche 
commenced (see Sec. II). Notice that the lifetime of an 
avalanche can range from seconds in the laboratory gran- 
ular piles to geological times in the landscape evolution. 
Therefore, the change of local stability conditions during 
the avalanche lifetime is a natural choice in the case of 
long relaxation times. A similar type of disorder in di- 
rected percolation processes was recently considered by 
Jensen |]T^ . 

We perform extensive numerical simulations for vari- 
ous values of the parameter /xg and lattice sizes L, and 
quantify the behavior by the landslide distributions of: 
(i) duration t — time that an instability lasts measured on 
the internal time scale; (ii) size s — area affected by an in- 
stability, and (iii) mass n — number of grains that exhibit 
slides during one avalanche, and (iv) by outflow current 
J — number of grains that fall off the open boundaries 
of the pile. Self-organized critical states are found for a 
range of values of the control parameter fio > ij-q ~ 0.4, 
which are characterized with multifractal scaling proper- 
ties and /iQ-dependent critical exponents. For fiQ < fi^ 
large discharging events occur occasionally, representing 
large-scale erosional reorganization of the system rather 
than fluctuations around a well defined critical state. 

The organization of the paper is as follows: In Sec. 

II we introduce the model and show two representative 
examples of landslides. The probability distributions of 
slides and their scaling properties are determined in Sec. 

III and IV for various values of the linear system size L 
and the parameter /io in the scaling region. Sec. V con- 
tains a short summary and the discussion of the results. 

II. MODEL AND LANDSLIDES 

We consider a square lattice oriented downward, with 
a dynamic variable, height h{i, j), associated to each site. 
The relaxation rules are a combination of (i) stochastic 
diffusion by two particles when h{i,j) > he with proba- 
bility which varies in time (see below); and (ii) de- 
terministic convection, when local slope <T{i,j) = h(i, j) — 
h(i + 1, j±) exceeds some critical value cr{i,j) > Uc- At 
each site the rule (ii) is applied by toppling one parti- 
cle along an unstable slope repeatedly until both local 
slopes drop below Uc- The system is updated in parallel. 



which leads to a well defined internal time scale of the 
relaxation process. The updating is stopped when all af- 
fected sites become temporarily stable. Here (i -I- 1, j±) 
are positions of two downward neighbors of the site(i, j). 
Mass flow is always downward, however, the instability 
can propagate backwards both due to nonlocal slope con- 
dition and due to time-dependent diffusion probability. 
We assume that diffusion probability fluctuates stochas- 
tically in time, but is space independent. Implementation 
of this rule is done as follows: We preset the threshold 
value /iQ which is the same for all sites in the system. 
Then at each site which is affected by an avalanche a 
new value ^,{t) is selected at each time step until the 
avalanche stops from set of random numbers evenly dis- 
tributed on the interval (0,1), and toppling is accepted 
if ^{t) < Ho, and rejected otherwise Therefore, for 
/io = 1 all sites topple (the rule becomes deterministic), 
whereas for /iq < 1 an unstable site might not topple at 
a given time t because of instantly low diffusion proba- 
bility p(t) = 1 — fi(t) < 1 — fiQ, however, it may topple 
at a later time step t' > t ii p{t') exceeds the threshold 
diffusion probability 1 — /xq . This temporally varying dis- 
order mimics changes in sticking properties with time, 
which then locally influence the angle of repose. This 
phenomenon can be of interest for the flow of granular 
materials with large effective friction, such as ricepiles 
in which the effects of granular boundaries may depend 
on the local dynamic variable h(i,j) and its derivatives. 
Therefore the difference /x(t) — /io is a measure of the dy- 
namic friction. Recently proposed models with stochastic 
critical slope rules in one dimension p5[ | proved very suc- 
cessful in describing the observed transport properties of 
ricepiles Whereas for avalanche distributions these 
models predict universal scaling exponents, in contrast 
to the experimental observations jsJI^-Q] . 

Another interesting example is represented by land- 
scape evolution, which can also be considered as a gran- 
ular flow iQl, in which local wetting properties fluctuate 
in time. By wetting, p{t) drops below the threshold dif- 
fusion probability 1 — /io, the grains stick together, and 
the system builds up large local slopes. At a later time 
t' these slopes may become unstable either when due to 
drying p{t') exceeds the threshold, or when the slopes 
become larger than critical. Two different classes of trig- 
gering mechanisms of landslides have been discussed in 
the literature : rainfall and water level which control 
soil moisture on one side, and ground motion, which leads 
to slope variation on the other. The values of measured 
exponents of landslide distributions are directly related 
to the locally prevailing triggering mechanism JTzt . In 
principle, threshold shear stress may depend on the slope 
angle and on soil properties, which are influenced by soil 
moisture. We assume that these two mechanisms are 
related dynamically. In the present model both mecha- 
nisms are effective: The soil moisture, which affects local 
height, varies stochastically in time at each site, whereas 
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we assume that the shear stress threshold depends only 
on the local angle and thus remains deterministic. More- 
over, by tuning the critical height mechanism via the pa- 
rameter /zq, we find nonuniversal critical properties and 
a transition to noncritical dynamic states, in qualitative 
agreement with experimental observations. A different 
model of landslides is obtained by "averaging out" the 
critical height mechanism and assuming stochastic vari- 
ations of critical slope, which can be viewed as one of 
few possible generalizations of stochastic critical slope 
models to two dimensions. So far the results of two- 
dimensional stochastic critical slope models are not avail- 
able in the literature [2^ . 

The system is perturbed by adding grains one at a time 
at a random site on the first row, thus increasing local 
height and slopes. Therefore, an instability (avalanche) 
can in principle start only from the top, however, sec- 
ondary avalanches are commencing from any affected site 
in the system, triggered either by high instant value of 
fi{t) or by supercritical slope. In order to have "clean" 
statistics, we start each avalanche from the top row and 
consider only those secondary avalanches which are spa- 
tially connected within certain correlation time tc- Here 
tc is not a prefixed parameter, but it is determined by 
the relaxation process itself. Typically tc is determined 
by the lifetime of the instability, thus ^ 1 for large 
relaxation events. There are two interesting limits of our 
model. In the limit /iq = 1 it reduces to the determin- 
istic directed model [Q, whereas for fj,o < 1 and in the 
limit when the correlation time is strictly equal to one, it 
reduces to the model considered in Ref. [ p^ . 

The temporally varying diffusion probability is a new 
ingredient of our model, which was not considered so far 
in CA models of granular flow. It appears to be responsi- 
ble both for new scaling properties and for the transition 
into the state dominated by large erosional avalanches. 
In Fig. 1 are shown two examples of simulated landslides 
with multiple topplings due to secondary avalanches up 
to fourth degree in the scaling region (top) and a large 
erosional event (bottom). 

III. PROBABILITY DISTRIBUTIONS OF SLIDES 
AND THEIR SCALING PROPERTIES 

In this section we present results of numerical simu- 
lations of avalanche statistics. As discussed in Sec. I, 
a landslide consists of many interpenetrating avalanches 
of different degree, which are spatially connected one to 
another within the life-time of the instability. For con- 
creteness, the probability distributions are determined 
for the whole relaxation event, which is equally termed 
as avalanche and/or landslide. We apply open boundary 
conditions in the perpendicular direction (see also later 
an example where periodic boundaries have been used). 
In most simulations we used he = 2 and (Tc = 8. By vary- 
ing the external parameter /uq between and 1 and lattice 



size L between 12 and 192, we determine the distributions 
of size, mass and duration of avalanches (slides). 

In Figs. 2 and 3 the distributions of avalanche dura- 
tion longer than t, P{t), size larger than s, -D(s), and 
mass larger than n, D{n), are shown for L = 128 and 
various values of the parameter /xq- (Notice that in the 
deterministic limit iiq — I the distributions D{s) and 
D{n) become identical, however, unbounded number of 
topplings at each site for fiQ < 1 leads to two distinct 
distributions.) For fj,Q < 1 a characteristic behavior 
with two scales appears: steep section corresponding to 
small avalanches, and flat section to large avalanches. 
The crossover length between small and large relaxation 
events varies with /uq, however, it remains small (cf. Figs. 
2 and 3), so that distributions of avalanches smaller than 
the crossover length extend only over one decade. Here 
we concentrate on the behavior of large avalanches (i.e., 
avalanches which arc larger than the crossover length). 
With lowering the threshold diffusion probability a 
large number of secondary instabilities develop, lead- 
ing to the flattening of the distributions. However, we 
find a power- law behavior P{t) i^"^*, D{s) - 3^-^% 
and D{n) ^ n^~'^^ , as long as /iq > 0.4. The expo- 
nents Tj, Ts and r„ appear to vary continuously with 
control parameter /ig, as shown in the insets to Figs. 
2 and 3. The character of the dynamics changes below 
/ig w 0.4, where only occasionally very large avalanches 
occur. We study in some more detail the relaxation clus- 
ters at /ig = 0.4. Numerical values of the exponents are 
Tt = 1.253, Ts = 1.202, and r„ = 1.190 for the distribu- 
tions of duration, size, and mass of avalanches, respec- 
tively. In addition, we have measured the distribution of 
linear elongation of avalanches in the direction of trans- 
port P{£) the mass-to-scale ratio with respect 
to parallel length {s)i ~ , and the average transverse 
extent {£±) - We find n = 1.578, D\\ = 1.572, and 
C = — 1 = 0.572 (estimated error bars ±0.03). These 
values are close to the numerical values of the exponents 
in the parity-conserving universality class |24| of branch- 
ing processes. On the other hand, the exponents govern- 
ing small events increase with decreasing /ig (cf. Fig. 2), 
reaching the values =1.92, tJ' — 1.67, and =1.45 
for the duration, size, and mass of small avalanches, re- 
spectively, at fiQ = /Iq. Notice that although the scale 
of the the distributions is small, being bounded by the 
crossover length, these values of the exponents indicate 
closeness of the mean-field universality class. 

A. Multifractal scaling properties of landslide 
distributions 

By varying the lattice size L with /ig fixed in the scaling 
region we study the finite-size effects on the distributions 
of avalanches. In contrast to the most of two-dimensional 
sandpile automata models in the literature, the present 
distributions do not obey simple finite-size scaling. In- 
stead, we find that different regions of a large avalanche 
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have different fractal properties and consequently their 
own exponents. The foUowing multifractal scahng form 

n 

P(X,L)^(L/Lo)^-("-) , (1) 

with 

„.(,„g|.)/(,og|.) (2) 

fits well our data with Lq — 1/4 and Xq = 1/4. (Here 
X = t,s, J). In Figs. 4 and 5 we show the probability 
distributions of duration and size, respectively, for five 
different lattice sizes L and for fixed fiQ = 0.7 . The cor- 
responding spectral functions (/it (at) vs. at and (/)s{as) 
vs. Qfs are shown in the insets to Figs. 4 and 5. 

IV. OUTFLOW CURRENT 

The outflow current results only from those avalanches 
that reach an open boundary of the system. The size of 
such events and their frequency is a relative measure of 
the transport processes which occur in the interior of the 
pile. The outflow current is easy to measure both in 
laboratory experiments and in natural landslides. For 
instance, the width of the sedimented layers of granular 
materials that occur below steep sections in mountains 
are directly related to the size of outflow avalanches from 
that section. Sensitivity of the outflow current distribu- 
tion P( J) to variations in the control parameter is mon- 
itored in our model for L — A8 with periodic boundary 
conditions in the perpendicular direction. In Fig. 6 we 
show the distribution P(J) vs. J for hq =1, 0.8, 0.6. , 
0.4 and 0.2 . Once again, the change in the character of 
the dynamics below (Iq is also seen in the outflow current, 
which becomes centered around a certain mean value (de- 
pending on the lattice size). Above /ig, we find that the 
outflow current distribution exhibits multifractal scaling 
properties according to Eqs. (0) and (H). The results 
for fiQ —0.7 and varying L from 12 to 192, obtained for 
open boundary conditions in perpendicular direction, are 
shown in Fig. 7. 

Additional information about transport processes in 
the interior of the system is obtained by measuring the 
outflow current as a function of time, and time intervals 
between successive outflow events. In the inset to Fig. 6 
we show the average time interval between outflow events 
as a function of the control parameter /io. The time in- 
tervals grow exponentially on lowering fiQ. In Fig. 8 the 
outflow current is shown as a function of time (measured 
on the external time scale, i.e., by the number of added 
particles), averaged over 1000 time steps for L =54 and 
with periodic perpendicular boundary conditions. For 
/io > 0.4 (cf. lower three panels), the outflow current 
fluctuates around mean value Jq = 1, thus balancing the 



input current and maintaining the steady states of the 
system (a steady state is characterized by balance be- 
tween input and output currents). The amplitude of the 
outflow events increases with decreasing /ig, and at the 
same time the frequency of events decreases. This be- 
havior is consistent with the histogram which is shown 
in Fig. 6. The character of the dynamics changes for 
/io < Mo (^^® top panel in Fig. 8), with dominating out- 
put events of large size and large time intervals between 
the events. At /zq = /io a dynamic phase transition oc- 
curs between critical steady states above /ig and states 
without long-range correlations below /io —Jip- (Similar 
phase transitions are found also in Refs. |l^ and [^ , 
however, in different universality classes.) Although for 
/io < /io the system is likely to build up a finite slope 
(unlimited piling is prevented by the deterministic criti- 
cal slope rule), preliminary results show that a substan- 
tial growth of the average slope occurs only for /io < 0.2, 
reaching the value Uc at jUq — > 0. Further work is neces- 
sary in order to investigate the universality class of this 
phase transition. 

V. DISCUSSION AND CONCLUSIONS 

In the present model, combined relaxation rules with 
temporal disorder are responsible for numerous after- 
avalanches, which lead to large relaxation events resem- 
bling sandslides in realistic granular materials. Numer- 
ical simulations show that such large relaxation events 
exhibit scaling behavior for a range of values of the the 
control parameter fJ-o ^ fJ-o ~ 0.4. The avalanche distri- 
butions are characterized by continuously varying scal- 
ing exponents, in qualitative agreement with the data 
collected from natural landslides. Moreover, comparison 
of the exponent of the avalanche mass distribution r„ 
for 0.4 < /io < 0.5 with the one that characterizes Hi- 
malayan sandslides reported in Ref. is satisfactory. 
For various lattice sizes the distributions are character- 
ized by multifractal rather than finite-size scaling proper- 
ties. The deterministic part of the relaxation rules leads 
to branching processes with, on the average, even number 
of offsprings. For this reason the scaling exponents for the 
distributions reach numerical values characteristic of the 
modulo-two conserving processes (also known as parity 
conserving processes) before scaling behavior disappears 
at fiQ = Mo- Below /ig the critical steady state is lost. The 
dynamics is dominated by large erosional avalanches in 
a region close to /ij and a net average slope appears for 
smaller values of fiQ. 
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FIG. 1. Two examples of avalanches running from left to 
right: (top) in the scaling region fio > Mo s^nd (below) in the 
region of erosional avalanches. Multiple topplings up to forth 
order are marked by different degrees of gray color. 
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FIG. 2. Double-logarithmic plot of the integrated prob- 
ability distribution of avalanche durations P{t) vs. duration 

t for L = 128 and various values of the control parameter 
/io=l, 0.9, 0.8, 0.7, 0.6, and 0.5 (top to bottom). Dashed and 
dotted lines indicate slopes of small and large avalanches, re- 
spectively. Inset: Scaling exponent of large avalanches n — 1 
vs. flo- 
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FIG. 3. Double-logarithmic plot of the integrated proba- 
bility distribution of avalanche size D{s) vs. size s (bottom) 
and mass D{n) vs. n (top), for L = 128 and for (top to bot- 
tom) /io=l, 0.8, 0.7, 0.6, and 0.5 . Inset: Scaling exponent of 
large avalanches Ts — 1 vs. fio (bottom figure) and Tn — 1 vs. 
/Uo (top figure). 
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FIG. 4. Double-logarithmic plot of the distribution P{t) 
vs. t for /uo=0.7 and for various lattice sizes L=12, 24, 48, 96, 
and 192 (left to right) with open boundary conditions. Inset: 
Multifractal scaling function (ptio-t) vs. at ■ 
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FIG. 5. Double-logarithmic plot of the distribution of sizes 
D(s) vs. s for the same set of parameters as Fig. 4. Inset: 
Multifractal scaling function 4's(as) vs. as ■ 



6 



8.0 



6.0 



a. 4.0 



2.0 



0.0 




0.0 1.0 2.0 3.0 



4.0 

In J 



5.0 6.0 7.0 



8.0 



FIG. 6. Double-logarithmic plot of the probability distri- 
bution of outflow current P{J) vs. J for L=48 with periodic 
boundary conditions in perpendicular direction, and for jiu = l, 
0.8, 0.6, 0.4, and 0.2 (top to bottom). Inset: Average time 
intervals between outflow events on the same lattice for crc=8 
(open triangles) and crc=4 (flUed triangles) . 
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FIG. 7. Distribution of outflow current measured with 
open boundary conditions for various lattice sizes L=12, 24, 
48, 96, and 192 (left to right) and for fixed //o=0.7. Inset: 
Spectral function (/>j(aj) vs. aj . 
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FIG. 8. Outflow current J(t) vs. time t (measured in 
number of added particles), averaged over 1000 time steps, 
for L=54 and Ato=l, 0.7, 0.4, and 0.3 (bottom to top) and 
with periodic boundary conditions. Dashed lines are mean 
values calculated by linear fits of the data for t > 40: (bot- 
tom to top) 0.9972, 1.0004, 0.9935 and 0.9892. Slopes of the 
dashed lines are smaller than 10~® in each case. 
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